﻿using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Linq;
using System.Text;
using Common;


namespace ThreeItems
{
    public class CentralizedAlg
    {
        /// <summary>
        /// 计算梯度时使用的算法
        /// </summary>
        public GradientType GradientMethod { get; set; }

        /// <summary>
        /// 计算步长使用的算法
        /// </summary>
        public StepLengthType StepLengthMethod { get; set; }

        /// <summary>
        /// 计算收益时使用的算法
        /// </summary>
        public ProfitType ProfitMethod { get; set; }

        /// <summary>
        /// 初始点类型
        /// </summary>
        public InitialPointType InitialPointMethod { get; set; }

        public double LinearSearchStep { get; set; }

        public double StepError { get; set; }

        public double Tolerences { get; set; }

        public double RepeatTimes { get; set; }

        //计算集中式模型
        public void RunComputing(ExampleInfo example)
        {
            double pnorm1 = (example.U1 + example.R12 * (example.U2 + example.O2) + example.R13 * (example.U3 + example.O3)) / (example.U1 + example.O1 + example.R12 * (example.U2 + example.O2) + example.R13 * (example.U3 + example.O3));
            example.MaxQC1 = NormDistHelper.QNorm(pnorm1, example.E1, example.SD1);
            double pnorm2 = (example.U2 + example.R21 * (example.U1 + example.O1) + example.R23 * (example.U3 + example.O3)) / (example.U2 + example.O2 + example.R21 * (example.U1 + example.O1) + example.R23 * (example.U3 + example.O3));
            example.MaxQC2 = NormDistHelper.QNorm(pnorm2, example.E2, example.SD2);
            double pnorm3 = (example.U3 + example.R31 * (example.U1 + example.O1) + example.R32 * (example.U2 + example.O2)) / (example.U3 + example.O3 + example.R31 * (example.U1 + example.O1) + example.R32 * (example.U2 + example.O2));
            example.MaxQC3 = NormDistHelper.QNorm(pnorm3, example.E3, example.SD3);

            //初始点 Q = 0
            if (InitialPointMethod == InitialPointType.Zero)
            {
                example.QC1 = 0;
                example.QC2 = 0;
                example.QC3 = 0;
            }
            else if (InitialPointMethod == InitialPointType.Classic)
            {
                example.QC1 = NormDistHelper.QNorm(example.UO1, example.E1, example.SD1);
                example.QC2 = NormDistHelper.QNorm(example.UO2, example.E2, example.SD2);
                example.QC3 = NormDistHelper.QNorm(example.UO3, example.E3, example.SD3);
            }

            //迭代
            bool stop = false;
            //存储连续小于容忍值的迭代
            int repeates = 0;

            double lastProfit = 0;

            for (int j = 0; j < 100; j++)
            {
                if (stop)
                {
                    break;
                }

                double profit = ComputeProfit(example);

                //todo:在此，记录每次迭代的数据，观察收敛效果
                IterateOnce(example);

                profit = ComputeProfit(example);

                double error = profit - lastProfit;
                lastProfit = profit;
                if (error < Tolerences)
                    repeates += 1;
                else //出现任意一次不能接受，重置计数器
                    repeates = 0;

                //判断是否连续N次可以接受
                if (repeates >= RepeatTimes)
                    stop = true;
            }

            double qc1 = example.QC1;
            double qc2 = example.QC2;
            double qc3 = example.QC3;

            Trace.WriteLine(string.Format("收益={0},QC1={1},QC2={2},QC3={3}", example.PC, qc1, qc2, qc3));

        }

        /// <summary>
        /// 一次迭代
        /// </summary>
        /// <param name="example"></param>
        /// <returns>误差e^2</returns>
        private void IterateOnce(ExampleInfo example)
        {
            //求梯度 二维数组
            double[] g = ComputeGradient(example);

            double qc1 = example.QC1;
            double qc2 = example.QC2;
            double qc3 = example.QC3;

            if (Moduel(g) == 0) //已经不能移动
            {
                return;
            }
            double[] ng = NormalizedGradient(g);


            //求步长
            double a = ComputeStepLength(example);
            //Trace.WriteLine("步长：" + a);
            //生成新的Q
            example.QC1 = qc1 + ng[0] * a;
            example.QC2 = qc2 + ng[1] * a;
            example.QC3 = qc3 + ng[2] * a;
        }

        #region 几个计算梯度的候选方法
        //计算梯度
        private double[] Gradient(ExampleInfo example)
        {
            //求梯度 二维数组
            var g = new double[3];

            double u1 = example.U1;
            double u2 = example.U2;
            double u3 = example.U3;

            double o1 = example.O1;
            double o2 = example.O2;
            double o3 = example.O3;

            double e1 = example.E1;
            double e2 = example.E2;
            double e3 = example.E3;

            double sd1 = example.SD1;
            double sd2 = example.SD2;
            double sd3 = example.SD3;

            double var1 = sd1 * sd1;
            double var2 = sd2 * sd2;
            double var3 = sd3 * sd3;

            double qc1 = example.QC1;
            double qc2 = example.QC2;
            double qc3 = example.QC3;


            double r12 = example.R12;
            double r13 = example.R13;
            double r21 = example.R21;
            double r23 = example.R23;
            double r31 = example.R31;
            double r32 = example.R32;


            //计算有效需求的期望和标准差
            double gama1 = (qc1 - ((qc1 - e1) * NormDistHelper.PNorm(qc1, e1, sd1) + var1 * NormDistHelper.DNorm(qc1, e1, sd1))) / e1;
            double gama2 = (qc2 - ((qc2 - e2) * NormDistHelper.PNorm(qc2, e2, sd2) + var2 * NormDistHelper.DNorm(qc2, e2, sd2))) / e2;
            double gama3 = (qc3 - ((qc3 - e3) * NormDistHelper.PNorm(qc3, e3, sd3) + var3 * NormDistHelper.DNorm(qc3, e3, sd3))) / e3;

            double e1_e = e1 - r21 * e2 * (1 - gama2) - r31 * e3 * (1 - gama3);
            double sd1_e = Math.Sqrt(var1 + Math.Pow(r21 * sd2 * (1 - gama2), 2) + Math.Pow(r31 * sd3 * (1 - gama3), 2));

            double e2_e = e2 - r12 * e1 * (1 - gama1) - r32 * e3 * (1 - gama3);
            double sd2_e = Math.Sqrt(var2 + Math.Pow(r12 * sd1 * (1 - gama1), 2) + Math.Pow(r32 * sd3 * (1 - gama3), 2));

            double e3_e = e3 - r13 * e1 * (1 - gama1) - r23 * e2 * (1 - gama2);
            double sd3_e = Math.Sqrt(var3 + Math.Pow(r13 * sd1 * (1 - gama1), 2) + Math.Pow(r23 * sd2 * (1 - gama2), 2));


            //需要确定相关系数矩阵
            var sigma12 = new double[2, 2]; // D_1 和D_2e的相关系数矩阵
            sigma12[0, 0] = var1;
            sigma12[0, 1] = sigma12[1, 0] = -r12 * (1 - gama1) * var1;
            sigma12[1, 1] = sd2_e * sd2_e;

            var sigma13 = new double[2, 2]; // D_1 和D_3e的相关系数矩阵
            sigma13[0, 0] = var1;
            sigma13[0, 1] = sigma13[1, 0] = -r13 * (1 - gama1) * var1;
            sigma13[1, 1] = sd3_e * sd3_e;

            var sigma21 = new double[2, 2]; //D_2和 D_1e 的相关系数矩阵
            sigma21[0, 0] = var2;
            sigma21[0, 1] = sigma21[1, 0] = -r21 * (1 - gama2) * var2;
            sigma21[1, 1] = sd1_e * sd1_e;

            var sigma23 = new double[2, 2]; //D_2和 D_3e 的相关系数矩阵
            sigma23[0, 0] = var2;
            sigma23[0, 1] = sigma23[1, 0] = -r23 * (1 - gama2) * var2;
            sigma23[1, 1] = sd3_e * sd3_e;

            var sigma31 = new double[2, 2]; //D_3和 D_1e 的相关系数矩阵
            sigma31[0, 0] = var3;
            sigma31[0, 1] = sigma31[1, 0] = -r31 * (1 - gama3) * var3;
            sigma31[1, 1] = sd1_e * sd1_e;

            var sigma32 = new double[2, 2]; //D_3和 D_2e 的相关系数矩阵
            sigma32[0, 0] = var3;
            sigma32[0, 1] = sigma32[1, 0] = -r32 * (1 - gama3) * var3;
            sigma32[1, 1] = sd2_e * sd2_e;

            //分别求两个梯度
            g[0] = r12 * (u2 + o2) *
                   NormDistHelper.PBvNorm(qc1, double.PositiveInfinity, double.NegativeInfinity, qc2, e1, e2_e, sigma12)
                   +
                   r13 * (u3 + o3) *
                   NormDistHelper.PBvNorm(qc1, double.PositiveInfinity, double.NegativeInfinity, qc3, e1, e3_e, sigma13)
                   + u1 - (u1 + o1) * NormDistHelper.PNorm(qc1, e1_e, sd1_e);
            g[1] = r21 * (u1 + o1) *
                   NormDistHelper.PBvNorm(qc2, double.PositiveInfinity, double.NegativeInfinity, qc1, e2, e1_e, sigma21)
                   +
                   r23 * (u3 + o3) *
                   NormDistHelper.PBvNorm(qc2, double.PositiveInfinity, double.NegativeInfinity, qc3, e2, e3_e, sigma23)
                   + u2 - (u2 + o2) * NormDistHelper.PNorm(qc2, e2_e, sd2_e);
            g[2] = r31 * (u1 + o1) *
                   NormDistHelper.PBvNorm(qc3, double.PositiveInfinity, double.NegativeInfinity, qc1, e3, e1_e, sigma31)
                   +
                   r32 * (u2 + o2) *
                   NormDistHelper.PBvNorm(qc3, double.PositiveInfinity, double.NegativeInfinity, qc2, e3, e2_e, sigma32)
                   + u3 - (u3 + o3) * NormDistHelper.PNorm(qc3, e3_e, sd3_e);


            return g;
        }

        //精确梯度计算 数值积分
        private double[] GradientNI(ExampleInfo example)
        {
            //求梯度 二维数组
            var g = new double[2];

            double u1 = example.U1;
            double u2 = example.U2;
            double o1 = example.O1;
            double o2 = example.O2;
            double e1 = example.E1;
            double e2 = example.E2;
            double sd1 = example.SD1;
            double sd2 = example.SD2;
            double var1 = sd1 * sd1;
            double var2 = sd2 * sd2;

            double qc1 = example.QC1;
            double qc2 = example.QC2;


            double r12 = example.R12;
            double r21 = example.R21;


            //计算有效需求的期望和标准差

            double e1_e = e1 - r21 * e2 + r21 * qc2;
            double sd1_e = Math.Sqrt(var1 + r21 * r21 * var2);

            double e2_e = e2 - r12 * e1 + r12 * qc1;
            double sd2_e = Math.Sqrt(var2 + r12 * r12 * var1);

            //需要确定相关系数矩阵
            //var corrs1 = new double[2, 2]; // D_1 和D_2e的相关系数矩阵
            //corrs1[0, 0] = 1;
            //corrs1[0, 1] = corrs1[1, 0] = -r12 * sd1 / sd2_e;
            //corrs1[1, 1] = 1;

            //var corrs2 = new double[2, 2]; //D_2和 D_1e 的相关系数矩阵
            //corrs2[0, 0] = 1;
            //corrs2[0, 1] = corrs2[1, 0] = -r21 * sd2 / sd1_e;
            //corrs2[1, 1] = 1;

            //联合概率密度
            double jointlyProb1 = (1 - NormDistHelper.PNorm(qc1, e1, sd1)) * NormDistHelper.PNorm(qc2, e2_e, sd1_e);

            //NormDistHelper.PBvNorm(qc1, double.PositiveInfinity, double.NegativeInfinity, qc2, e1,
            //                                         e2_e, corrs1);
            double jointlyProb2 = (1 - NormDistHelper.PNorm(qc2, e2, sd2)) * NormDistHelper.PNorm(qc1, e1_e, sd2_e);
            //NormDistHelper.PBvNorm(qc2, double.PositiveInfinity, double.NegativeInfinity, qc1, e2,
            //                                         e1_e, corrs2);

            //singleProb 使用全概率公式来计算，也就是计算两个条件概率
            double singleProb1 = NormDistHelper.PNorm(qc1, e1, sd1) * NormDistHelper.PNorm(qc2, e2, sd2)
                                 + (1 - NormDistHelper.PNorm(qc2, e2, sd2)) * NormDistHelper.PNorm(qc1, e1_e, sd1_e);
            // jointlyProb2;
            double singleProb2 = NormDistHelper.PNorm(qc1, e1, sd1) * NormDistHelper.PNorm(qc2, e2, sd2)
                + (1 - NormDistHelper.PNorm(qc1, e1, sd1)) * NormDistHelper.PNorm(qc2, e2_e, sd2_e);
            //jointlyProb1;

            //分别求两个梯度
            g[0] = r12 * (u2 + o2) * jointlyProb1 + u1 - (u1 + o1) * singleProb1;
            g[1] = r21 * (u1 + o1) * jointlyProb2 + u2 - (u2 + o2) * singleProb2;


            return g;
        }

        public double[] NormalizedGradient(double[] g)
        {
            double module = Math.Sqrt(Math.Pow(g[0], 2) + Math.Pow(g[1], 2) + Math.Pow(g[2], 2));
            if (module != 0)
            {
                return new[] { g[0] / module, g[1] / module, g[2] / module };
            }
            else throw new ArgumentException("模不能是0");
        }

        ////计算梯度 归一化的
        //public double[] NormalizedGradient(ExampleInfo example)
        //{
        //    return NormalizedGradient(Gradient(example));
        //}

        //求梯度的模
        public double Moduel(double[] g)
        {
            return Math.Sqrt(Math.Pow(g[0], 2) + Math.Pow(g[1], 2) + Math.Pow(g[2], 2));
        }

        #endregion


        #region 几个计算步长的候选方法

        //简单方法 固定增长 直到找到收益开始下降
        private double StepLengthSimple(ExampleInfo example)
        {
            double qc1 = example.QC1;
            double qc2 = example.QC2;

            double[] g = ComputeGradient(example);
            double[] ng = NormalizedGradient(g);

            double a0;
            double lastProfit = ComputeProfit(example);

            //先顺序找出初始步长
            for (int i = 1; ; i++)
            {
                example.QC1 += LinearSearchStep * ng[0];
                example.QC2 += LinearSearchStep * ng[1];

                //夹角法计算导数
                double profit = ComputeProfit(example);

                if (profit <= lastProfit)
                {
                    a0 = LinearSearchStep * (i - 0.5); //和上次步长取平均数
                    // Debug.WriteLine(string.Format("找到初始步长： {0}", a0));
                    break;
                }
                else
                {
                    lastProfit = profit;
                }
            }
            //}
            //Trace.WriteLine("初始步长：" + a0);

            //恢复QC
            example.QC1 = qc1;
            example.QC2 = qc2;

            return a0;
        }

        //使用Wolfe-Powell条件进行不精确一维搜索
        private double StepLengthPowell(ExampleInfo example)
        {
            //计算梯度
            double[] g = ComputeGradient(example);
            double[] d = NormalizedGradient(g);

            double a = 0.1;
            double b = 0.7;
            double step = 10;

            double qc1 = example.QC1;
            double qc2 = example.QC2;

            double profit = ComputeProfit(example);
            //寻找步长
            while (true)
            {

                example.QC1 = qc1 + a * step * d[0];
                example.QC2 = qc2 + a * step * d[1];

                if (ComputeProfit(example) + 0.1 >= (profit + a * step * (g[0] * d[0] + g[1] * d[1])))//满足1
                {
                    double[] g2 = ComputeGradient(example);
                    if (g2[0] * d[0] + g2[1] * d[1] <= b * (g[0] * d[0] + g[1] * d[1]))//满足2
                    {
                        break;
                    }
                    else
                    {
                        step = step * 3 / 2;
                    }
                }
                else
                {
                    step = step / 2;
                }
            }

            //恢复QC
            example.QC1 = qc1;
            example.QC2 = qc2;

            return step;
        }

        //0.618法
        private double StepLengthGoldenSection(ExampleInfo example)
        {
            double qc1 = example.QC1;
            double qc2 = example.QC2;
            double qc3 = example.QC3;

            //计算梯度
            double[] g = ComputeGradient(example);
            double[] d = NormalizedGradient(g);

            //double a1 = random.Next();
            //double b1 = random.Next();
            double a = 0;// Math.Max(a1, b1);
            double b = 1000;// Math.Min(a1, b1);
            double err = 0.0001;

            double t = (Math.Sqrt(5) - 1) / 2; // = 0.618

            double u = a + (1 - t) * (b - a);
            double v = a + t * (b - a);

            while (true)
            {
                if ((b - a) < err)
                    break;
                else
                {
                    example.QC1 = qc1 + d[0] * u;
                    example.QC2 = qc2 + d[1] * u;
                    example.QC3 = qc3 + d[2] * u;
                    double profitLeft = ComputeProfit(example);

                    example.QC1 = qc1 + d[0] * v;
                    example.QC2 = qc2 + d[1] * v;
                    example.QC3 = qc3 + d[2] * v;
                    double profitRight = ComputeProfit(example);

                    if (profitLeft < profitRight)
                    {
                        a = u;
                        u = v;
                        v = a + t * (b - a);
                    }
                    else
                    {
                        b = v;
                        v = u;
                        u = a + (1 - t) * (b - a);
                    }
                }
            }

            //恢复QC
            example.QC1 = qc1;
            example.QC2 = qc2;
            example.QC3 = qc3;

            return (a + b) / 2;
        }

        //牛顿迭代法 - 计算步长
        private double StepLengthNewton(ExampleInfo example)
        {
            double a;
            /*
             寻找步长之前已经有一个确定的方向g
             在这个方向上使用 Newton–Raphson (牛顿-拉夫逊) 方法求解一阶导数为0这个方程，得到的解就是步长
             
             */

            double u1 = example.U1;
            double u2 = example.U2;
            double o1 = example.O1;
            double o2 = example.O2;
            double e1 = example.E1;
            double e2 = example.E2;
            double sd1 = example.SD1;
            double sd2 = example.SD2;
            double var1 = sd1 * sd1;
            double var2 = sd2 * sd2;

            double qc1 = example.QC1;
            double qc2 = example.QC2;


            double r12 = example.R12;
            double r21 = example.R21;


            //计算有效需求的期望和标准差
            double pnorm1 = NormDistHelper.PNorm(qc1, e1, sd1);
            double gama1 = pnorm1 +
                           (qc1 * (1 - pnorm1) - (var1) * NormDistHelper.DNorm(qc1, e1, sd1)) /
                           e1;

            double pnorm2 = NormDistHelper.PNorm(qc2, e2, sd2);
            double gama2 = pnorm2 +
                           (qc2 * (1 - pnorm2) - var2 * NormDistHelper.DNorm(qc2, e2, sd2)) /
                           e2;

            double e1_e = e1 - r21 * e2 * (1 - gama2);
            double sd1_e = Math.Sqrt(var1 + Math.Pow(r21 * sd2 * (1 - gama2), 2));

            double e2_e = e2 - r12 * e1 * (1 - gama1);
            double sd2_e = Math.Sqrt(var2 + Math.Pow(r12 * sd1 * (1 - gama1), 2));


            double[] g = ComputeGradient(example);
            double[] ng = NormalizedGradient(g);
            // List<double> steps = new List<double>();//存储找到的步长

            #region 测试 观察收益对于步长的导数的变化

            /*
            double step = 0.1;
            int points = 90;
            double[] dev1s = new double[points];
            double[] dev2s = new double[points];
            double[] steps = new double[points];
            double[] profits = new double[points];
            for (int j = 0; j < points; j++)
            {
                double u1_t = u1;
                double u2_t = u2;
                double o1_t = o1;
                double o2_t = o2;



                //double dev_1 = ng[0] * (u1_t + o1_t) * NormDistHelper.PNorm(qc1 + step * ng[0], e1_e, sd1_e)
                //        + ng[1] * (u2_t + o2_t) * NormDistHelper.PNorm(qc2 + step * ng[1], e2_e, sd2_e)
                //        - u1_t * ng[0]
                //        - u2_t * ng[1];

                ////double dev_2 = ng[0] * NormDistHelper.DNorm(qc1 + step * ng[0], e1_e, sd1_e)
                ////           + ng[1] * NormDistHelper.DNorm(qc2 + step * ng[1], e2_e, sd2_e);

                //double dev_2 = ng[0] * ng[0] * (u1 + o1) * NormDistHelper.DNorm(qc1 + step * ng[0], e1_e, sd1_e)
                //       + ng[1] * ng[1] * (u2 + o2) * NormDistHelper.DNorm(qc2 + step * ng[1], e2_e, sd2_e);
                steps[j] = step;
                
                 profits[j] = Profit(example);

                 example.QC1 += step * ng[0];
                 example.QC2 += step * ng[1];
                 step += 0.1;
                //dev1s[j] = dev_1;
                //dev2s[j] = dev_2;
            }
            //RHelper.DrawPlot(steps, dev1s);
            //RHelper.DrawPlot(steps, dev2s);
            RHelper.DrawPlot(steps, profits);
             */

            #endregion

            double a0 = 0;
            //double lastDevirant1 = ng[0] * (u1 + o1) * NormDistHelper.PNorm(qc1, e1_e, sd1_e)
            //                      + ng[1] * (u2 + o2) * NormDistHelper.PNorm(qc2, e2_e, sd2_e)
            //                      - u1 * ng[0]
            //                      - u2 * ng[1];

            double lastDevirant = ng[0] * g[0] + ng[1] * g[1];


            //在此输出夹角法和老师推导的 看看之间差距有多大？如果不大就继续用这个方法。。。。


            //if (lastDevirant < 0)//todo:有没有可能步长为0时导数是大于0的，也就是需要向梯度反方向移动。要不就是有可能要不就是需要检查程序错误
            //{

            var exam = new ExampleInfo();
            //先顺序找出初始步长
            for (int i = 1; ; i++)
            {

                exam.QC1 = qc1 + LinearSearchStep * i * ng[0];
                exam.QC2 = qc2 + LinearSearchStep * i * ng[1];

                exam.SD1 = example.SD1;
                exam.SD2 = example.SD2;

                exam.E1 = example.E1;
                exam.E2 = example.E2;

                exam.U1 = example.U1;
                exam.U2 = example.U2;

                exam.O1 = example.O1;
                exam.O2 = example.O2;

                exam.R21 = example.R21;
                exam.R12 = example.R12;

                double[] g2 = ComputeGradient(exam);


                //计算方程左侧的值(lambda导数)
                //double devirant1 = ng[0] * (u1 + o1) * NormDistHelper.PNorm(qc1 + linearSearchStep * i * ng[0], e1_e, sd1_e)
                //                  + ng[1] * (u2 + o2) * NormDistHelper.PNorm(qc2 + linearSearchStep * i * ng[1], e2_e, sd2_e)
                //                  - u1 * ng[0]
                //                  - u2 * ng[1];


                //夹角法计算导数
                double devirant = ng[0] * g2[0] + ng[1] * g2[1];

                //Trace.WriteLine(string.Format("{0}\t{1}", devirant, -devirant1));

                //注意下面的e1_e和sd1_e，需要修正
                if (devirant == 0)
                {
                    a0 = LinearSearchStep * i;
                    //Debug.WriteLine(string.Format("找到初始步长： {0}", a0));
                    break;
                }

                if (devirant * lastDevirant <= 0)
                {
                    a0 = LinearSearchStep * (i - 0.5); //和上次步长取平均数
                    // Debug.WriteLine(string.Format("找到初始步长： {0}", a0));
                    break;
                }
                else
                {
                    lastDevirant = devirant;
                    //Debug.WriteLine(string.Format("第{0}次寻找初始步长，a0={1}", i, linearSearchStep * i));
                }
            }
            //}
            //Trace.WriteLine("初始步长：" + a0);

            return a0;

            //x_k出的梯度 和 x_k+1出的梯度的数量积 是导数
            //这个牛拉方法 要求初始点很严格
            //计算初始点
            //double a1_1 = (NormDistHelper.QNorm(u1 / (u1 + o1), e1_e, sd1_e) - qc1) / ng[0];
            //double a2_1 = (NormDistHelper.QNorm(u2 / (u2 + o2), e2_e, sd2_e) - qc2) / ng[1];

            //a1 = Math.Min(a1_1 ,a2_1);
            if (a0 == 0) //todo:check 比0稍大一点的初始步长
            {
                a = 0;
            }
            else
            {
                double a1 = a0;
                double a2; //a1 上一个 a2当前,也是最优步长
                for (int i = 0; ; i++)
                {
                    //double[] g2 = Gradient(example);
                    // a2 = a1 - Profit(example) / Math.Sqrt((g[0] * g2[0] + g[1] * g2[1]));//在梯度方向的导数值是梯度的模


                    //函数值 fx
                    double fx = ng[0] * (u1 + o1) * NormDistHelper.PNorm(qc1 + a1 * ng[0], e1_e, sd1_e)
                                + ng[1] * (u2 + o2) * NormDistHelper.PNorm(qc2 + a1 * ng[1], e2_e, sd2_e)
                                - u1 * ng[0]
                                - u2 * ng[1];


                    //导数 fxd
                    double fxd = ng[0] * ng[0] * (u1 + o1) * NormDistHelper.DNorm(qc1 + a1 * ng[0], e1_e, sd1_e)
                                 + ng[1] * ng[1] * (u2 + o2) * NormDistHelper.DNorm(qc2 + a1 * ng[1], e2_e, sd2_e);

                    //todo:test this : double fxd = ng[0] * NormDistHelper.DNorm(qc1 + a1 * ng[0], e1_e, sd1_e)
                    //            + ng[1] * NormDistHelper.DNorm(qc2 + a1 * ng[1], e2_e, sd2_e);
                    //推导公式
                    a2 = a1 - fx / fxd;
                    //Debug.WriteLine(string.Format("fx={0},fxd={1},x2={2}", fx, fxd, a2));

                    //steps.Add(a2);
                    //Debug.WriteLine(string.Format("第{0}次寻找步长，为{1}", i + 1, a2));

                    //比较
                    if (Math.Abs(a1 - a2) < StepError)
                    {
                        a = a2;
                        break;
                    }
                    else
                    {
                        //example.QC1 += (a2 - a1) * ng[0];
                        //example.QC2 += (a2 - a1) * ng[1];
                        a1 = a2;
                    }
                }
            }
            Trace.WriteLine("牛拉步长：" + a);

            /**/
            return a;
        }

        #endregion


        #region 目标收益计算

        //计算利润
        private double Profit(ExampleInfo example)
        {
            double u1 = example.U1;
            double u2 = example.U2;
            double u3 = example.U3;

            double o1 = example.O1;
            double o2 = example.O2;
            double o3 = example.O3;

            double qc1 = example.QC1;
            double qc2 = example.QC2;
            double qc3 = example.QC3;



            //计算有效需求的期望和标准差
            ComputeApproximateDemand(example);

            double e1_e = example.E1_E;
            double sd1_e = example.SD1_E;

            double e2_e = example.E2_E;
            double sd2_e = example.SD2_E;

            double e3_e = example.E3_E;
            double sd3_e = example.SD3_E;


            double profit1 = u1 * qc1 -
                             (u1 + o1) * ((qc1 - e1_e) * NormDistHelper.PNorm(qc1, e1_e, sd1_e) +
                              sd1_e * sd1_e * NormDistHelper.DNorm(qc1, e1_e, sd1_e));


            double profit2 = u2 * qc2 -
                             (u2 + o2) *
                             ((qc2 - e2_e) * NormDistHelper.PNorm(qc2, e2_e, sd2_e) +
                              sd2_e * sd2_e * NormDistHelper.DNorm(qc2, e2_e, sd2_e));

            double profit3 = u3 * qc3 -
                             (u3 + o3) *
                             ((qc3 - e3_e) * NormDistHelper.PNorm(qc3, e3_e, sd3_e) +
                              sd3_e * sd3_e * NormDistHelper.DNorm(qc3, e3_e, sd3_e));

            double profit = profit1 + profit2 + profit3;

            example.PC1 = profit1;
            example.PC2 = profit2;
            example.PC3 = profit3;
            example.PC = profit;

            return profit;
        }

        //计算利润 数值积分
        private double ProfitNI(ExampleInfo example)
        {

            double u1 = example.U1;
            double u2 = example.U2;
            double o1 = example.O1;
            double o2 = example.O2;
            double e1 = example.E1;
            double e2 = example.E2;
            double sd1 = example.SD1;
            double sd2 = example.SD2;

            double qc1 = example.QC1;
            double qc2 = example.QC2;

            double r12 = example.R12;
            double r21 = example.R21;

            double var1 = sd1 * sd1;
            double var2 = sd2 * sd2;


            double e1_e = e1 - r21 * e2 + r21 * qc2;
            double sd1_e = Math.Sqrt(var1 + r21 * r21 * var2);
            double var1_e = sd1_e * sd1_e;

            double e2_e = e2 - r12 * e1 + r12 * qc1;
            double sd2_e = Math.Sqrt(var2 + r12 * r12 * var1);
            double var2_e = sd2_e * sd2_e;

            ////singleProb 使用全概率公式来计算，也就是计算两个条件概率
            //double singleProb1 = NormDistHelper.PNorm(qc2, e2, sd2) * NormDistHelper.PNorm(qc1, e1, sd1)
            //                     + (1 - NormDistHelper.PNorm(qc2, e2, sd2)) * NormDistHelper.PNorm(qc1, e1_e, sd1_e);
            //// jointlyProb2;
            //double singleProb2 = NormDistHelper.PNorm(qc1, e1, sd1) * NormDistHelper.PNorm(qc2, e2, sd2)
            //    + (1 - NormDistHelper.PNorm(qc1, e1, sd1)) * NormDistHelper.PNorm(qc2, e2_e, sd2_e);


            double p1 = NormDistHelper.PNorm(qc1, e1, sd1);
            double p2 = NormDistHelper.PNorm(qc2, e2, sd2);
            double p1_e = NormDistHelper.PNorm(qc1, e1_e, sd1_e);
            double p2_e = NormDistHelper.PNorm(qc2, e2_e, sd2_e);

            double d1 = NormDistHelper.DNorm(qc1, e1, sd1);
            double d2 = NormDistHelper.DNorm(qc2, e2, sd2);
            double d1_e = NormDistHelper.DNorm(qc1, e1_e, sd1_e);
            double d2_e = NormDistHelper.DNorm(qc2, e2_e, sd2_e);

            //D_ie < Q_i
            double prob1 = p2 * p1 + (1 - p2) * p1_e;
            double prob2 = p1 * p2 + (1 - p1) * p2_e;

            //期望值计算 D_ie < Q_i
            double prob1_e = p2 * (e1 * p1 - var1 * d1) + (1 - p2) * (e1_e * p1_e - var1_e * d1_e);
            double prob2_e = p1 * (e2 * p2 - var2 * d2) + (1 - p1) * (e2_e * p2_e - var2_e * d2_e);

            //double prob1_e = e1_e * p1_e - var1_e * d1_e;
            //double prob2_e = e2_e * p2_e - var2_e * d2_e;

            double expect1 = qc1 * prob1 - prob1_e;
            double expect2 = qc2 * prob2 - prob2_e;


            //double expect1 = qc1 * (p2 * p1 + (1 - p2) * p1_e) - p2 * (e1 * p1 - var1 * d1) - (1 - p2) * (e1_e * p1_e - var1_e * d1_e);
            //double expect2 = qc2 * (p1 * p2 + (1 - p1) * p2_e) - p1 * (e2 * p2 - var2 * d2) - (1 - p1) * (e2_e * p2_e - var2_e * d2_e);

            //在这里 分开求 试试看。。
            //。。。。。


            //分别求Profit1 和 Profit2
            double profit1 = u1 * qc1 - (u1 + o1) * expect1;
            double profit2 = u2 * qc2 - (u2 + o2) * expect2;

            //NormDistHelper.PNorm(qc2, e2, sd2)*

            double profit = profit1 + profit2;

            example.PC1 = profit1;
            example.PC2 = profit2;
            example.PC = profit;

            return profit;
        }

        //计算近似需求 使用结构来传递数据
        public void ComputeApproximateDemand(ExampleInfo example)
        {
            double e1 = example.E1;
            double e2 = example.E2;
            double e3 = example.E3;

            double sd1 = example.SD1;
            double sd2 = example.SD2;
            double sd3 = example.SD3;

            double var1 = sd1 * sd1;
            double var2 = sd2 * sd2;
            double var3 = sd3 * sd3;

            double qc1 = example.QC1;
            double qc2 = example.QC2;
            double qc3 = example.QC3;


            double r12 = example.R12;
            double r13 = example.R13;
            double r21 = example.R21;
            double r23 = example.R23;
            double r31 = example.R31;
            double r32 = example.R32;


            //计算有效需求的期望和标准差
            double gama1 = (qc1 - ((qc1 - e1) * NormDistHelper.PNorm(qc1, e1, sd1) + var1 * NormDistHelper.DNorm(qc1, e1, sd1))) / e1;
            double gama2 = (qc2 - ((qc2 - e2) * NormDistHelper.PNorm(qc2, e2, sd2) + var2 * NormDistHelper.DNorm(qc2, e2, sd2))) / e2;
            double gama3 = (qc3 - ((qc3 - e3) * NormDistHelper.PNorm(qc3, e3, sd3) + var3 * NormDistHelper.DNorm(qc3, e3, sd3))) / e3;

            example.E1_E = e1 - r21 * e2 * (1 - gama2) - r31 * e3 * (1 - gama3);
            example.SD1_E = Math.Sqrt(var1 + Math.Pow(r21 * sd2 * (1 - gama2), 2) + Math.Pow(r31 * sd3 * (1 - gama3), 2));

            example.E2_E = e2 - r12 * e1 * (1 - gama1) - r32 * e3 * (1 - gama3);
            example.SD2_E = Math.Sqrt(var2 + Math.Pow(r12 * sd1 * (1 - gama1), 2) + Math.Pow(r32 * sd3 * (1 - gama3), 2));

            example.E3_E = e3 - r13 * e1 * (1 - gama1) - r23 * e2 * (1 - gama2);
            example.SD3_E = Math.Sqrt(var3 + Math.Pow(r13 * sd1 * (1 - gama1), 2) + Math.Pow(r23 * sd2 * (1 - gama2), 2));
        }

        #endregion

        //计算梯度 选择器
        public double[] ComputeGradient(ExampleInfo example)
        {
            if (GradientMethod == GradientType.Pricise)
                return GradientNI(example);
            if (GradientMethod == GradientType.Approximation)
                return Gradient(example);
            throw new ArgumentException("传入的梯度计算方法参数不正确");
        }

        //计算步长 选择器
        public double ComputeStepLength(ExampleInfo example)
        {
            if (StepLengthMethod == StepLengthType.GoldenSection)
                return StepLengthGoldenSection(example);
            if (StepLengthMethod == StepLengthType.PowellCondition)
                return StepLengthPowell(example);
            throw new ArgumentException("传入的梯度计算方法参数不正确");
        }

        //计算收益 选择器
        public double ComputeProfit(ExampleInfo example)
        {
            if (ProfitMethod == ProfitType.Pricise)
                return ProfitNI(example);
            if (ProfitMethod == ProfitType.Approximation)
                return Profit(example);
            throw new ArgumentException("传入的梯度计算方法参数不正确");
        }
    }

    /// <summary>
    /// 初始点类型
    /// </summary>
    public enum InitialPointType
    {
        //0为初始点
        Zero = 1,

        //经典多产品报童模型的最优解为初始点 交叉系数为0
        Classic = 2
    }

    public enum GradientType
    {
        //近似梯度
        Approximation = 1,

        //精确梯度
        Pricise = 2
    }

    public enum ProfitType
    {
        //近似梯度
        Approximation = 1,

        //精确梯度
        Pricise = 2
    }

    public enum StepLengthType
    {
        //牛顿拉弗逊方法 
        NewtonRaphson = 1,

        //黄金分割发
        GoldenSection = 2,

        //Powell条件 非精确搜索
        PowellCondition = 3,
    }
}
